Algorytm automatycznej hp adaptacji został po raz pierwszy szczegółowo opisany w książkach prof. Leszka Demkowicza, matematyka polskiego pochodzenia, pracującego na Uniwersytecie Teksańskim w Austin. Analizę właściwości matematycznych algorytmu hp-adaptacji dokonał prof. Ivo Babuška, amerykanin czeskiego pochodzenia, również pracujący na Uniwersytecie Teksańskim w Austin [1] [2] [3] [4].

Algorytm hp-adaptacji zilustrowany jest na rysunkach Rys. 1 - Rys. 5.
Konwencja przyjęta podczas rysowania siatek obliczeniowych jest następująca. Kolory na rysunkach Rys. 1, Rys. 3 i Rys. 5 oznaczają stopnie wielomianów rozpiętych na krawędziach i wnętrzach elementów. Na każdej krawędzi ustalony jest jeden stopień wielomianu, natomiast na każdym wnętrzu elementu ustalone są stopień w kierunku poziomym oraz stopień w kierunku pionowym.
Na przykład, na rysunku Rys. 1 wszystkie wielomiany na wszystkich krawędziach mają stopień 2, oraz wszystkie wielomiany na wnętrzach wszystkich elementów mają stopień 2 w kierunku pionowym i poziomym. Na rysunku Rys. 3 wszystkie wielomiany na wszystkich krawędziach mają stopień 3, oraz wszystkie wielomiany na wnętrzach wszystkich elementów mają stopień 3 w kierunku pionowym i poziomym. Z koleji na rysunku Rys. 5 przeglądając elementy wierszami od lewej do prawej strony, w pierwszym wierszu na pierwszym elemencie wszystkie wielomiany na krawędziach mają stopień 2, z wyjątkiem wielomianu na dolnej krawędzi elementu, który ma stopień 1. Na pierwszym elemencie wielomiany we wnętrzu mają stopień 2 w kierunku pionowym i poziomym. Na drugim elemencie, na krawędziach lewej, górnej i dolnej, wielomiany mają stopień 2, natomiast we wnętrzu stopień 3 w kierunku pionowym i stopień 2 w kierunku poziomym. Na trzecim i czwartym elemencie wszystkie wielomiany na krawędziach mają stopień 3. Podobnie wielomiany we wnętrzach w kierunkach pionowym i poziomym. W drugim wierszu elementów, pierwszy element (czyli piąty globalnie) posiada wielomiany stopnia 3 na lewej i prawej krawędzi, wielomian stopnia 1 na górnej i dolnej krawędzi. We wnętrzu element ten posiada wielomiany stopnia 1 w kierunku poziomym i wielomiany stopnia 3 w kierunku pionowym. Drugi element w drugim wierszu posiada stopień wielomianu 2 w kierunku poziomym oraz stopień wielomianu 3 w kierunku pionowym. Na trzecim i czwartym elemncie wszystkie wielomiany na wszystkich krawędziach mają stopień 3, oraz wszystkie wielomiany na wnętrzach wszystkich elementów mają stopień 3 w kierunku pionowym i poziomym. Pozostałe elementy w trzecim i czwartym wierszu mają wielomiany rozłożone symetrycznie w stosunku do wielomianów w pierwszym i drugim wierszu na pierwszym i drugim elemencie.
W jaki sposób dokonywany jest wybór optymalnego sposobu adaptacji pojedynczego elementu? Zauważmy że w przypadku algorytmu hp adaptacji mamy wiele możliwości modyfikacji pojedynczego elementu:

1. Pozostawienie elementu bez zmian
2. Możliwe jest połamanie elementu na dwa nowe elementy w kierunku pionowym
3. Możliwe jest połamanie elementu na dwa nowe elementy w kierunku poziomym
4. Możliwe jest połamanie elementu na cztery nowe elementy (dwa w kierunku poziomym i dwa w kierunku pionowym)

Ponadto dla każdego połamanego elementu możliwa jest modyfikacja stopnia wielomianów we wnętrzu elementu

1. Pozostawienie stopni elementu bez zmian
2. Możliwe jest podniesienie stopnia we wnętrzu elementu w kierunku poziomym \( (p_h,p_v)\rightarrow (p_h+1,p_v) \)
3. Możliwe jest podniesienie stopnia we wnętrzu elementu w kierunku pionowym \( (p_h,p_v)\rightarrow (p_h,p_v+1) \)
4. Możliwe jest podniesienie stopnia we wnętrzu elementu w kierunku poziomym i pionowym \( (p_h,p_v)\rightarrow (p_h+1,p_v+1) \)

Stopnie na krawędziach zmodyfikowanych elementów ustalane są na podstawie reguły minimum. Innymi słowy stopień wielomianu na krawędzi dzielonej przez dwa sąsiadujące elementy ustalany jest jako minimum stopni we wnętrzu elementów w odpowiednim kierunku.
Dla krawędzi pionowej otoczonej przez dwa elementy \( (p^1_h,p^2_v) \) oraz \( (p^2_h,p^2_v) \) ustalamy stopień \( p=\textrm{min} \{ p^1_v,p^2,v\} \).
Dla krawędzi poziomej otoczonej przez dwa elementy \( (p^1_h,p^2_v) \) oraz \( (p^2_h,p^2_v) \) ustalamy stopień \( p=\textrm{min} \{ p^1_h,p^2,h\} \).
Mamy więc bardzo dużo możliwości adaptacji pojedynczego elementu (szacując na podstawie wymienionych możliwych rodzajów adaptacji mamy 4 (modyfikacje stopnia elementu bez łamania elementu) + 4*4 (połamanie elementu na dwa elementy w kierunku poziomym i po cztery możliwości modyfikacji stopnia każdego elementu) + 4*4 (połamanie elementu na dwa elementy w kierunku pionowym i po cztery możliwości modyfikacji stopnia każdego elementu) + 4*4*4*4 (połamanie elementu na cztery elementy i po cztery możliwości modyfikacji stopnia każdego elementu). W sumie 4+4*4+4*4+4*4*4*4=292 możliwości.
W jaki sposób podejmowana jest decyzja o wyborze rodzaju adaptacji pojedynczego elementu?
Decyzja podejmowana jest zgodnie z następującym algorytmem.

Innymi słowy, wybieramy taki rodzaj adaptacji dla elementu, który pozwala uzyskać największy spadek błędu najmniejszym kosztem. Wielkość ta reprezentowana jest przez prędkość spadku błędu, która rośnie wraz ze spadkiem błędu, ale maleje wraz z poniesionym kosztym (z liczbą funkcji dodaną do elemetu, ponieważ koszt obliczenia na danym elemencie zależy od liczby niewiadomych, które są współczynnikami funkcji bazowych). Algorytm ten zilustrowany jest na rysunku Rys. 6

Dla przykładowego problemu transportu ciepła na obszarze w kształcie litery L możliwe jest uzyskanie siatki hp adaptacyjnej która pozwala na rozwiązanie problemu z błędem względnym 0.0001. Siatka taka przedstawiona jest na rysunkach Rys. 7 - Rys. 13

Algorytm automatycznej hp adaptacji jest więc bardzo skomplikowany i trudny do zaimplementowania. Dlaczego więc wart jest on naszego zainteresowania? Rozważmy możliwe różne algorytmy adaptacyjne dla przykładowego problemu transportu ciepła na obszarze w kształcie litery L, czyli

1. Algorytm jednorodnej h adaptacji, w którym równomiernie na całym obszarze zwiększana jest liczba elementów
2. Algorytm jednorodnej p adaptacji, używający siatki zbudowanej z 3 elementów i zwiększający równomiernie stopień wielomianów we wnętrzach wszystkich elementów w kierunku poziomym oraz pionowym, oraz stopień wielomianów na wszystkich krawędziach
3. Algorytm jednorodnej p adaptacji, używający siatki zbudowanej z 12 elementów i zwiększający równomiernie stopień wielomianów we wnętrzach wszystkich elementów w kierunku poziomym oraz pionowym, oraz stopień wielomianów na wszystkich krawędziach
4. Algorytm automatycznej p adaptacji, zwiększający stopień wielomianów we wnętrzach wybranych elementów w wybranych kierunkach, i modyfikujący stopnie na krawędziach zgodnie z regułą minimum
5. Algorytm automatycznej h adaptacji, łamiący wybrane elementy w wybranych kierunkach
6. Algorytm automatycznej hp adaptacji, opisany w tym rozdziale

Jeśli narysujemy wykres zbieżności tych algorytmów w taki sposób, że na osi poziomej umieścimy rozmiar siatki rozumiany jako liczbę funkcji bazowych na całej siatce, co odpowiada liczbie niewiadomych współczynników tych funkcji bazowych czyli rozmiarowi macierzy układu równań do rozwiązania, a na osi pionowej błąd wzgłędny liczony w normie H1 pomiędzy rozwiązaniem na siatce rzadkiej i siatce gęstej
\( \|u_{h,p}-u_{h/2,p+1}\|_{H^1}=\frac{\int_{\Omega} (u_{h,p}-u_{h/2,p+1})^2+(\frac{\partial (u_{h,p}-u_{h/2,p+1}}{\partial x})^2+(\frac{\partial (u_{h,p}-u_{h/2,p+1}}{\partial x})^2 }{\int_{\Omega} (u_{h/2,p+1})^2+(\frac{\partial (u_{h/2,p+1}}{\partial x})^2+(\frac{\partial (u_{h/2,p+1}}{\partial x})^2 } \)
wówczas okaże się że jedynie algorytm automatycznej hp adaptacji posiada cechę zbieżności eksponencjalnej. Jest on zdecydowanie najszybszy, i potrafi rozwiązać zadany problem obliczeniowy z dokładnością niemożliwą do uzyskania przez inne algorytmy adaptacyjne.

W praktyce większość problemów inżynieryjnych wymaga dokładności z zakresu 5 procent i wielomianów kwadratowych. Są jednak takie zadania obliczeniowe, gdzie wysoka dokładność rozwiązania jest niezbędna. Przykładem takich zadań obliczeniowych są na przykład symulacje przepływów wokół skrzydeł samolotów, które wymagają bardzo bardzo dużej dokładności w warstwie przyściennej na skrzydłach samolotu, lub symulacje propagacji fal elektromagnetzcznych w warstwach górotworu, w celu identyfikacji złóż ropy naftowej, które wymagają bardzo dużej dokładności w okolicy anteny odbiorczej.